## Visualization of the results
# newdat <- data.frame(
# AGE = rep(18:90, length.out = 500),
# female = rep(0:1, length.out = 500),
# black = rep(0:1, length.out = 500),
# hispanic = rep(0:1, length.out = 500),
# own = rep(0:1, length.out = 500),
# pid = rep(-1:1, length.out = 500),
# ideo = rep(-1:1, length.out = 500),
# college = rep(0:1, length.out = 500),
# married = rep(0:1, length.out = 500),
# income_gt50k = rep(0:1, length.out = 500),
# region = rep(c("Northeast", "South", "Midwest", "West"), length.out = 500))
########## Simulation for need-based policy #################
newdat <- subset(df_model, select = -c(ideo)) 
newdat <- cbind(newdat, predict(m1, newdat, type = "probs"))

lnewdat <- melt(newdat, id.vars = c("AGE", "female", "pid"),
                variable.name = "Level", value.name="Probability")

library(tidyr)
lnewdat <- pivot_longer(newdat, cols = 11:15)
names(lnewdat)[11] <- "Support" 
names(lnewdat)[12] <- "PredictedProb"

lnewdat$pid <- car::recode(as.numeric(lnewdat$pid),
                           "-2='All';
                                        -1='Democrats';
                                        0='Independents';
                                        1='Republicans'")

lnewdat$income_gt50k <- car::recode(as.numeric(lnewdat$income_gt50k),
                                    "0='Below Median Income level';
                                        1='Above Median Income level'")

lnewdat$black <- car::recode(as.numeric(lnewdat$black),
                             "0='Not being an African american';
                                        1='Being an African american'")

lnewdat<-lnewdat[!is.na(lnewdat$pid)& !is.na(lnewdat$income_gt50k) & !is.na(lnewdat$black),]

# lnewdat$Support <- car::recode(as.numeric(lnewdat$Support),
#                                    "0='Strongly oppose';
#                                    .25='Somewhat oppose';
#                                    .5='Neither support nor oppose';
#                                    .75='Somewhat support';
#                                    1='Strongly support'")
# 
# lnewdat$Support <- as.factor(lnewdat$Support) 
# lnewdat$Support <- factor(lnewdat$Support,
#                               levels = c("Strongly support", 
#                                          "Somewhat support",
#                                          "Neither support nor oppose",
#                                          "Somewhat oppose",
#                                          "Strongly oppose"))


lnewdat_2 <- lnewdat %>%
  group_by(Support, black, pid, income_gt50k) %>% 
  summarise(PredictedProbMean = mean(PredictedProb, na.rm = TRUE))
lnewdat_2$Support <- as.numeric(lnewdat_2$Support) 


lnewdat_2 %>%
  # filter(black == 0 & hispanic == 0 & region == "Northeast" & own == 0 & college == 0 & married == 0 & female == 0) %>%
  # filter(pid == NA) %>%
  ggplot(., aes(x = Support, y = PredictedProbMean, color = factor(black))) +
  geom_line() + 
  facet_grid(income_gt50k ~ pid) + 
  ggtitle("Predicted Probability for Different Support Levels \n by Political Party, Income, and Ethnicity ") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(color = "Ethnicity")+
  xlab("Support Level") + 
  ylab("Predicted probability of a person at a given support level")

#################################################################################### 
#################### redistributive treatment funding ################################ 
#################################################################################### 

newdat <- subset(df_model, select = -c(ideo)) 
newdat <- cbind(newdat, predict(m_ability2, newdat, type = "probs"))

lnewdat <- melt(newdat, id.vars = c("AGE", "female", "pid"),
                variable.name = "Level", value.name="Probability")

library(tidyr)
lnewdat <- pivot_longer(newdat, cols = 11:15)
names(lnewdat)[11] <- "Support" 
names(lnewdat)[12] <- "PredictedProb"

lnewdat$pid <- car::recode(as.numeric(lnewdat$pid),
                           "-2='All';
                                        -1='Democrats';
                                        0='Independents';
                                        1='Republicans'")

lnewdat$income_gt50k <- car::recode(as.numeric(lnewdat$income_gt50k),
                                    "0='Below Median Income level';
                                        1='Above Median Income level'")

lnewdat$black <- car::recode(as.numeric(lnewdat$black),
                             "0='Not being an African american';
                                        1='Being an African american'")

lnewdat<-lnewdat[!is.na(lnewdat$pid)& !is.na(lnewdat$income_gt50k) & !is.na(lnewdat$black),]

# lnewdat$Support <- car::recode(as.numeric(lnewdat$Support),
#                                    "0='Strongly oppose';
#                                    .25='Somewhat oppose';
#                                    .5='Neither support nor oppose';
#                                    .75='Somewhat support';
#                                    1='Strongly support'")
# 
# lnewdat$Support <- as.factor(lnewdat$Support) 
# lnewdat$Support <- factor(lnewdat$Support,
#                               levels = c("Strongly support", 
#                                          "Somewhat support",
#                                          "Neither support nor oppose",
#                                          "Somewhat oppose",
#                                          "Strongly oppose"))


lnewdat_2 <- lnewdat %>%
  group_by(Support, black, pid, income_gt50k) %>% 
  summarise(PredictedProbMean = mean(PredictedProb, na.rm = TRUE))
lnewdat_2$Support <- as.numeric(lnewdat_2$Support) 


lnewdat_2 %>%
  # filter(black == 0 & hispanic == 0 & region == "Northeast" & own == 0 & college == 0 & married == 0 & female == 0) %>%
  # filter(pid == NA) %>%
  ggplot(., aes(x = Support, y = PredictedProbMean, color = factor(black))) +
  geom_line() + 
  facet_grid(income_gt50k ~ pid) + 
  ggtitle("Predicted Probability for Different Support Levels \n by Political Party, Income, and Ethnicity ") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(color = "Ethnicity")+
  xlab("Support Level") + 
  ylab("Predicted probability of a person at a given support level")

#################################################################################### 
#################### building of a treatment facility nearby (distance_support) ###### 
#################################################################################### 

newdat <- subset(df_model, select = -c(ideo)) 
newdat <- cbind(newdat, predict(m_dist2, newdat, type = "probs"))

lnewdat <- melt(newdat, id.vars = c("AGE", "female", "pid"),
                variable.name = "Level", value.name="Probability")

library(tidyr)
lnewdat <- pivot_longer(newdat, cols = 11:15)
names(lnewdat)[11] <- "Support" 
names(lnewdat)[12] <- "PredictedProb"

lnewdat$pid <- car::recode(as.numeric(lnewdat$pid),
                           "-2='All';
                                        -1='Democrats';
                                        0='Independents';
                                        1='Republicans'")

lnewdat$income_gt50k <- car::recode(as.numeric(lnewdat$income_gt50k),
                                    "0='Below Median Income level';
                                        1='Above Median Income level'")

lnewdat$black <- car::recode(as.numeric(lnewdat$black),
                             "0='Not being an African american';
                                        1='Being an African american'")

lnewdat<-lnewdat[!is.na(lnewdat$pid)& !is.na(lnewdat$income_gt50k) & !is.na(lnewdat$black),]

# lnewdat$Support <- car::recode(as.numeric(lnewdat$Support),
#                                    "0='Strongly oppose';
#                                    .25='Somewhat oppose';
#                                    .5='Neither support nor oppose';
#                                    .75='Somewhat support';
#                                    1='Strongly support'")
# 
# lnewdat$Support <- as.factor(lnewdat$Support) 
# lnewdat$Support <- factor(lnewdat$Support,
#                               levels = c("Strongly support", 
#                                          "Somewhat support",
#                                          "Neither support nor oppose",
#                                          "Somewhat oppose",
#                                          "Strongly oppose"))


lnewdat_2 <- lnewdat %>%
  group_by(Support, black, pid, income_gt50k) %>% 
  summarise(PredictedProbMean = mean(PredictedProb, na.rm = TRUE))
lnewdat_2$Support <- as.numeric(lnewdat_2$Support) 


lnewdat_2 %>%
  # filter(black == 0 & hispanic == 0 & region == "Northeast" & own == 0 & college == 0 & married == 0 & female == 0) %>%
  # filter(pid == NA) %>%
  ggplot(., aes(x = Support, y = PredictedProbMean, color = factor(black))) +
  geom_line() + 
  facet_grid(income_gt50k ~ pid) + 
  ggtitle("Predicted Probability for Different Support Levels \n by Political Party, Income, and Ethnicity ") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(color = "Ethnicity")+
  xlab("Support Level") + 
  ylab("Predicted probability of a person at a given support level")